Distributions of forces and 'hydrodynamic clustering' in a shear thickened colloid. 
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Sheared concentrated colloids with short range polymer coats are examined via simulations. Dis- 
tributions of force are found to be sums of exponentials. The 'hydrodynamic clustering' underlying 
the shear thickening effect is shown, in this system, to be a network of percolating coat contacts, 
with coats both in compression and extension. The geometry and kinetics of this network are ex- 
plored along with its relation to the bulk stress tensor. Particles experience strong fluctuations in 
force over epochs circa 10% strain. Density fluctuations up to glassy volume fractions give rise to 
scattering peaks at low Q in the thickening regime. There is some commonality in the physics of 
this colloid particulate system and other granular media. 

Pacs Numbers: 83.50.-v, 47.20.Ft, 83.50.Qm, 83.70.Hq 



What is the physics of concentrated particles in flow 
? Detailed observation is hard, however, simulations of 
powders/granular media |^ have provided much in- 
sight. They have shown exponential distributions of force 
and force-chains. The origin of the distribution of force 
has motivated much theory [3-6], although in the main 
for static systems. 

This letter examines force transmission in steady states 
of sheared colloid systems. It concerns shear thickening, 
a dramatic effect in which, at high shear rates, concen- 
trated colloids develop high, and sometimes diverging, 
viscosities [Q. Understanding thickening is of consider- 
able interest [8-12]. The concept of 'hydrodynamic clus- 
tering' has been introduced to explain thickening, but its 
precise meaning has not been explored. The main aim 
of this letter is to clarify 'hydrodynamic clustering' , but 
on route to this commonalities will be found with dry 
powders. 

This work considers colloids of rigid particles sus- 
pended in a fluid. Particle motion is overdamped. The 
equation of motion is 
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where V,Fc;Fb are 6iV velocity/angular- velocity and 
colloid and Brownian forces/torque vectors and R is a 
6A^ X QN drag matrix of dissipative hydrodynamic inter- 
actions. 

The results below are, in the main, for a model which 
includes just the squeeze mode of hydrodynamic lubrica- 
tion forces between near neighbours. For spheres i and 
j with velocities Vi and vj this determines a force fy on 
particle i 



fij = -a(/i)(vi - Vj) -nyny, 



(2) 



where ny is the unit vector from centre i to j, a{h) the 
drag coefficient and h is the gap between the sphere sur- 
faces. R is made up of such pair terms. In the absence of 
dissipations for shear modes, the particles are 'friction- 
less'. Some simulations have also been performed that 
include the shear mode. 



The algorithm for simulating Eq. (1) with pair terms 
is given in [|l^ . The lubrication force opposes the relative 
motion of the particles. At the high concentrations of in- 
terest it is argued |l^ that these terms will dominate the 
long ranged mobility matrix in shear flow. In the case of 
a hard sphere model the drag coefficient a(h) is given by 
the well known Reynolds formula which diverges as 1 /h. 
The particles simulated here have crude models for short 
range polymer coats with both a conservative repulsive 
force and a lubrication interaction modified from that of 
hard spheres. The coats define a particle contact which 
is crucial to the physics discussed below. 

Experiment [^,and simulation [|l^, [|l^ establish that 
thickening is dominated by contributions from the lubri- 
cation forces. However, Reynolds lubrication is insufh- 
cient to drive strong thickening. Crude models of particle 
pairs shearing past each other jisj , lead to the conclusion 
that any divergence of bulk properties can only be as the 
logarithm of the inverse gap of closest approach hmin ■ It 
was suggested [^ that N-body effects may be necessary 
to explain the stronger divergences of experiment. 

However, simulations of spheres with Reynolds lubrica- 
tion [|lj] [Q do show only logarithmic divergence, much 
weaker than experiment [|ll| . The correct interpretation 
of the arguments in rcf. |15| is not that N-body effects are 
necessary, but that lubrication interactions stronger than 
Reynolds are necessary. Lubrication interaction can be 
strongly enhanced by polymer coats on the particles [ p^ . 
On compression of the coats the lubrication divergence 
switches to a higher power of (l//i) than Reynolds lu- 
brication. The strength of its lubrication coefficient, de- 
pends on its porosity (the smaller the pores, the stronger 
the lubrication coefficient; see for a plot of this in- 
teraction). The coat also has an elastic modulus. It was 
shown previously ||l^ that this model fits experimental 
data jn). 

Simple pair arguments [ pT| explain the onset of thick- 
ening. Thickening occurs once the bulk stress is sufficient 
to compress the coats against their conservative (spring) 
forces to a point at which the relaxation time of the coats. 
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Tc is longer than the inverse shear rate I/7. This is con- 
firmed by the simulations. 

Although we can understand the onset of thickening at 
the pair level, there is also a dependence on volume frac- 
tion, with strong thickening occurring only for (j}y > 0.50. 
Simulations do reveal many body effects. Simulations 
motivated a model of thickening by hydrodynamic 
clustering: agglomerates of particles bound together by 
lubrication forces. Assuming such clusters are rigid, ex- 
plains enhanced stress. However lubrication forces are 
only active under relative motion, so rigidity and agglom- 
eration under lubrication are not obviously compatible. 
Indeed the 'clusters' are reported to be transient [ 18| and 
observed to rapidly disappear after cessation of flow ^ - 
they exist in response to the stress. 

In the simulations below, the thickening regime when 
7Tc > 1 is identified with a network of coat contacts of 
mean coordination in the range 5-6. It is this network 
that constitutes 'hydrodynamic clustering'. The geome- 
try and kinetics of this network will be explored in detail. 

The simulated model has a coat thickness set at O.Old 
and other parameters choosen to fit experiments on per- 
spex particles [0. The shear rate is defined as the di- 
mensionless Peclet number. Figure 1 shows simulation 
data at volume fraction 0^, = 0.54. Below Pe = 200, 
the viscosity is close to flat and the particles flow in an 
ordered array of strings aligned along the flow direction. 
The system jumps from ordered flow to disordered flow 
from Pe = 200 to Pe — 500 where it continues to be dis- 
ordered and thicken at higher rates. This order-disorder 
transition during the initial stages of thickening is com- 
mon for mono-disperse spheres, however it is not suffi- 
cient nor necessary for thickening which is a feature of the 
disordered flow. We report elsewhere [ p^ systems where 
the onset of thickening is from regimes of disorder and 
charged systems at (f)y < 0.50 which show order disorder 
changes with minimal thickening. In any case, this let- 
ter is a study of the steady state thickening effect within 
the disordered regime. (The steady state should be con- 
trasted with jamming in hard spheres po|] . A jamming 
at a critical shear stress can also occur for models with 
coat interactions if the springs have a maximal force) 

Thickening is accompanied by large normal stress dif- 
ferences. In agreement with experiment, A^i is negative 
pl| and the coat lubrication forces are the dominant con- 
tribution to the stress tensor. At Pe = 200, jTc > 1 and 
remains so up to circa Pe = 10^ where the spring force 
of the highly compressed coats grows rapidly and the di- 
vergence of the lubrication coefficient weakens, bringing 
jTc close to unity and lowering the thickening. 

Figure 2 shows, for several shear rates, the probabil- 
ity distribution for the magnitude of the lubrication force 
on individual bonds. Bonds are defined as Voronoi near 
neighbours. Data is shown at four different shear rates. 
It was confirmed that the effects reported are indepen- 
dent of systems size above N=500 and independent of 



time step over a decade of variation. On the linear-log 
plot the distributions at high force are clearly seen to be 
exponential. 
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FIG. 1. Viscosity vs dimensionless shear rate (Peclet). 
Simulation at (f>v — 0.54 of colloid spheres with model 
polymer coats of thickness 0.01 diameters and porisity 
dp = O.OOObdiam.. Open symbols squeeze lubrication only, 
closed symbols with shear terms. Circles A'^ = 50 particles, 
squares A'^ = 200 and triangles A'' = 2000. Inset normal stress 
differences Ni/Pe (open circles) and N2/Pe (open diamonds). 

In the ordered phase at Pe = 100 the distribution 
is a single exponential for the bonds in extension and 
dominated by a single exponential for the bonds in com- 
pression (there is negligible 2nd exponential in the tail). 
These are not observable on the scale of figure 2. 

In the thickening regime at Pe = 500, 1000, 3000 the 
distributions over the full range are fit by the sum of 
two exponential decays, one particular example is shown 
in the inset. The system is characterised by two well 
separated characteristic forces. Bulk averages and the 
thickening are dominated by the growth of the high force 
distributions. Examination shows that the distribution 
with the decay at higher force (henceforth the contact 
network) is just that of the bonds with coats in con- 
tact; 45-50% of the bonds are in the contact network; 
the mean number of contacts per particles is 5.2, 5.7, 5.9 
for Pe = 500, 1000, 3000 respectively It is noted that the 
coordination number is approaching, from below, that of 
an isostatic network of frictionless particles In the 
ordered regime a percolating contact network also exists, 
but is of lower coordination circa 3. 

The finding of an exponential distribution of forces is 
common to recent simulations [Qand experiment p^ ] 
on dry granular media. In shaken powders a bimodal dis- 
tribution was found In sheared powders large fluc- 
tuations in normal stress were reported ||2^. Note, we 
have also investigated force distributions when aggrega- 
tion forces are present and for hard spheres with 
Reynolds lubrication and short repulsive springs. Both 
these show exponential distributions. The distribution 
is not the result of a particular force law. Physical ar- 
gument for this distribution have been given [p[; it is 
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geometric in origin. 
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FIG. 2. Probability density of forces (units kT/d) on 
bonds. From the most inner pair to the most outer pair 
Pe = 500, 1000 and 3000 . Inset, a log-log plot of the dis- 
tribution of force magnitudes for the compressive bonds at 
Pe — 1000 (e.g. the circles in the main plot 2nd in from the 
left). Also shown in the inset is the fit of the whole distribu- 
tion too a sum of two exponentials (solid line) and the fit of 
just the high force part (dashed line) to a single exponential. 



The geometry of the network was studied. Overall, 
the network is branched. By thresholding on force or 
interparticle gap, clusters of M particles were defined 
and their radius of gyration tensor examined. At high 
force above the high characteristic forces, these had one 
eigenvalue scaling as M and others scaling roughly root 
M. Although linear at high force, they are not straight 
'force-chains' Q. By summing just over the contacts par- 
tial stress tensors, a, and a partial fabric tensor, F, were 
studied. 



F 



(3) 



For the systems at Pe=500,1000 and 3000, table 1 
gives data from the high force distributions. For the 
Pc=100 system it shows data for bonds with \fbond\ > 
2 X lO'^. The flow geometry has y(x) the shear gra- 
dient(flow) direction with (x = —y) the compression 
and (x = y) the extension directions. The table shows 
the viscosity Uxy/Pe and first normal stress difference 
-^1 = {(^xx — o'yy), and the x and y components of the 
normalised principal eigenvector of the fabric tensor. The 
ordered flow at Pe=100 has contacts in compression and 
extension just close to the gradient direction. For the 
thickened states contacts both in compression and ex- 
tension contribute to the shear stress. However, contacts 
in extension/compression contribute roughly 2/3 to 1/3 
of the normal stress difference. For contacts under com- 
pression, the fabric tensor has its principal axis just below 
the compression direction and is relatively independent 



of shear rate. For the bonds in extension its principal axis 
lies above the extension direction, thus enhancing their 
contribution to iVi, it tends more to the extension direc- 
tion the higher the shear rate. The lubrication forces are 
symmetric from extension to compression, but the spring 
forces are not. It is the latter force that is responsible 
for the asymmetry between the fabric of the compression 
and extension wings. 
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Table 1 Contents are described in the text. 
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FIG. 3. Short history of the most negative (compression) 
and most positive (extension) eigenvalues of the stress tensor 
summed just for the bonds surrounding one particular particle 
in the simulations at Pe=3000. The epoch shown is from 
strain 50 to 52. 



The kinetics of the networks was examined. It is 
straightforward to define an effective deformation ten- 
sor for a cluster. 'Rigid' clusters of particles defined on 
force or gap thresholds would have very low extension 
and compression (relative to the applied shear) and just 
rotate in the flow. No clear evidence of these was found, 
although some large clusters of contacts did have eigen- 
values for the symmetric part of their deformation tensors 
as low as 1/4 of that for affine shear. Figure 3 shows the 
large fluctuations in stress on a particle during flow; they 
build up and decay over epochs of 10% or so strain, at 
the applied shear rate this is much faster than the natural 
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relaxation of the coat. They are driven by the flow that 
is changing the geometrically determined stress concen- 
tration in the system. (In the frame co-rotating with the 
flow the systems experience a changing biaxial stress.) 
Large changes in the distribution of forces under a 'small 
change' in the direction of applied stress has been termed 
'fragility' It is unclear if the 10% strain represents a 
'small change' and whether the changes in the network 
are an example of 'fragility'. 

A key question is whether the contact networks have 
an inherent length scale. The data does suggest some 
correlation with density fluctuations. The structure fac- 
tor in the thickening regime has peaks oriented along the 
flow direction inside the near neighbour ring. Similar 
peaks were observed in Neutron scattering from thick- 
ened colloids § . Selecting particles with high forces and 
computing a partial structure factor also resulted in low 
Q-peaks oriented along the flow direction. Typically a 
particle must have candidate bonds for high forces close 
to the shear plane. If these are to carry high forces, and 
if the bonds through the particle are not co-linear, the 
other bonds must be able to supply large lateral forces to 
maintain overall balance. A locally high density in which 
particles coats were compressed in all directions would 
provide this. The distribution of a local volume fraction 
parameter defined as a sphere volume divide by the vol- 
ume of its Voronoi cell was examined. Over all particles 
this is peaked at 0.54. If the distribution is computed just 
over particles selected on a high threshold of force magni- 
tude, the peak is shifted to around local volume fractions 
0.58-0.6, circa the hard sphere glass transition. Only to 
this degree can the network be considered 'clustering'. 
However, some particles at low local volume fractions 
can also carry high forces. 

This work has detailed the nature of steady state thick- 
ening in colloids. Some features were common with the 
physics of slow powder flow - despite the different na- 
ture of the physical interactions. Clearly this is due to 
the generic importance of contact geometry in particle 
flows. Hopefully this will lead to exchanges of ideas and 
results. It is noted that the force network here does not 
conform to the simple notion of 'force chains' prevalent 
in much of the literature, it is branched and coordinated 
close to the iso-static limit, but very high force parts were 
like short segments of random walks. The existence and 
relevance of exponential distributions of force have been 
questioned, but figure 2 shows that, at least in flow, this 
is a significant feature. 

For the colloid community, the results have given a 
much needed clarification and detailing of the nature of 
'hydrodynamic clustering'. There are no rigid clusters 
as some have assumed, but a network of contacts whose 
natural relaxation is slower than the shear time. This 
network flows with large fluctuations in lubrication force 
driven by the changing geometry of contacts in the shear 
flow. It had not been realised that bonds both in com- 



pression and extension are involved. The origin of normal 
stress differences is due to the fabric of the contact net- 
work. Low Q density fluctuations observed in experiment 
are also associated with the contact network. The growth 
of the high force distribution and the percolating contact 
network could be taken as an 'order-parameter' for the 
thickening transition. 
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